#include "PSI1D_Recycling.h"
#include "../Parallel/ParallelMPI.h"
#include "../Field/Field2D.h"
#include "../Species/Species1D.h"

#include <cmath>
#include <iomanip>
#include <algorithm>
#include <ostream>
#include <sstream>

using namespace std;


//Constructor
PSI1D_Recycling::PSI1D_Recycling(PicParams* params_, ParallelMPI* pmpi_, vector<Species*>* vecSpecies_, PSIParameter psi_param_, int psi_number_):
    PSI1D(params_, pmpi_, vecSpecies_, psi_param_, psi_number_)
{
    log_info<<"psi type: "<<psi_param.psi_type<<"  psi position: "<<psi_param.psi_position;

}

PSI1D_Recycling::~PSI1D_Recycling()
{

}


void PSI1D_Recycling::perform_psi(int itime)
{
    //ion
    auto s1 = static_cast<Species1D*>((*vecSpecies)[psi_param.species_number_list[0]]);

    //atom
    auto s2 = static_cast<Species1D*>((*vecSpecies)[psi_param.species_number_list[1]]);


    auto& psi_particle1 = s1->psi_particles_;

    //cout<<"recycliong begein "<<endl;

    if(psi_param.psi_position == "left" && pmpi->isWestern())
    {
        //cout<<"recycliong begein left"<<endl;
        n_particle_emit_left_ = 0;
        for(auto& particle: psi_particle1)
        {
            if(particle.position[0] < 0.0)
            {
                n_particle_emit_left_++;
            }
        }

        double n_particle_emit_left_double = n_particle_emit_left_ * psi_param.recycling_factor;
        n_particle_emit_left_ = n_particle_emit_left_double;
        n_particle_emit_left_rem_ += (n_particle_emit_left_double - n_particle_emit_left_);
        if(n_particle_emit_left_rem_ >= 1.0)
        {
            n_particle_emit_left_++;
            n_particle_emit_left_rem_ -= 1.0;
        }

        auto cell_particles = s2->particles_.front();
        int i_particle_begin = cell_particles->size();
        cell_particles->add_particles(n_particle_emit_left_);
        //if(n_particle_emit_left_ > 0) cout<<"recycliong begein left "<<n_particle_emit_left_<<endl;
        for(int i_particle = i_particle_begin; i_particle < cell_particles->size(); i_particle++)
        {
            auto& particle = cell_particles->data_[i_particle];
            particle.position[0] = (((double)rand() / RAND_MAX)) * params->cell_length[0] * psi_param.position_offset;

            double ran;
            do {
                ran = (double)rand() / RAND_MAX;
            }
            while (ran == 0.0);
            //Velocity magnitude: from Maxwell velocity distribution
            //The angle between velocity of emitted particle and the surface normal: cosine
            //     cos(alpha) = sqrt(random number 0-1)
            //The azimuthal angle is uniformly distributed on the interval [0 2pi]
            double psm = sqrt(2.0 * params->const_e * psi_param.recycling_temperature / s2->species_param.mass) * sqrt(-log(ran));
            double cosAlpha = sqrt((double)rand() / RAND_MAX);
            double sinAlpha = sqrt(1.0 - cosAlpha * cosAlpha);
            double phi   = 2.0 * M_PI*(double)rand() / RAND_MAX;

            particle.velocity[0] = abs(psm * cosAlpha);
            particle.velocity[1] = psm * sinAlpha * cos(phi);
            particle.velocity[2] = psm * sinAlpha * sin(phi);
        }

    }
    else if(psi_param.psi_position == "right" && pmpi->isEastern())
    {
        //cout<<"recycliong begein right "<<psi_particle1.size()<<endl;
        n_particle_emit_right_ = 0;
        for(auto& particle: psi_particle1)
        {
            //cout<<"recycliong begein right000"<<particle.position[0]<<" "<<params->sim_length[0]<<endl;
            if(particle.position[0] > params->sim_length[0])
            {
                n_particle_emit_right_++;
            }
        }
        //if(n_particle_emit_right_ > 0) cout<<"recycliong begein right111 "<<n_particle_emit_right_<<endl;

        double n_particle_emit_right_double = n_particle_emit_right_ * psi_param.recycling_factor;
        n_particle_emit_right_ = n_particle_emit_right_double;
        n_particle_emit_right_rem_ += (n_particle_emit_right_double - n_particle_emit_right_);
        if(n_particle_emit_right_rem_ >= 1.0)
        {
            n_particle_emit_right_++;
            n_particle_emit_right_rem_ -= 1.0;
        }

        auto cell_particles = s2->particles_.back();
        int i_particle_begin = cell_particles->size();
        cell_particles->add_particles(n_particle_emit_right_);
        //if(n_particle_emit_right_ > 0) cout<<"recycliong begein right "<<n_particle_emit_right_<<endl;
        for(int i_particle = i_particle_begin; i_particle < cell_particles->size(); i_particle++)
        {
            auto& particle = cell_particles->data_[i_particle];
            particle.position[0] = params->sim_length[0] - (((double)rand() / RAND_MAX)) * params->cell_length[0] * psi_param.position_offset;

            double ran;
            do {
                ran = (double)rand() / RAND_MAX;
            }
            while (ran == 0.0);
            //Velocity magnitude: from Maxwell velocity distribution
            //The angle between velocity of emitted particle and the surface normal: cosine
            //     cos(alpha) = sqrt(random number 0-1)
            //The azimuthal angle is uniformly distributed on the interval [0 2pi]
            double psm = sqrt(2.0 * params->const_e * psi_param.recycling_temperature / s2->species_param.mass) * sqrt(-log(ran));
            double cosAlpha = sqrt((double)rand() / RAND_MAX);
            double sinAlpha = sqrt(1.0 - cosAlpha * cosAlpha);
            double phi   = 2.0 * M_PI*(double)rand() / RAND_MAX;

            particle.velocity[0] = -abs(psm * cosAlpha);
            particle.velocity[1] = psm * sinAlpha * cos(phi);
            particle.velocity[2] = psm * sinAlpha * sin(phi);
        }
    }
    else
    {
        
    }

}

